import pandas as pd
import openpyxl
import math
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.mlab as mlab
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
from scipy.optimize import curve_fit

# Constants
N = 5  # Number of data points to average
# Read data from Excel file
title = "measurement_data"
df = pd.read_excel(title + '.xlsx')
title = "Sample_1"
# Extract data from DataFrame
mass = df['m'].tolist()
time = df['t'].tolist()

# Parameters
mass_MSS_start=df.iloc[0, 2]
mass_MSS_end=df.iloc[1, 2]
D=df.iloc[2,2]
area=df.iloc[3, 2]
T_dry=df.iloc[4,2]
Psi=df.iloc[5,2]
air_flow1=df.iloc[6, 2]
air_flow2=df.iloc[7, 2]
T_wet=df.iloc[8,2]
air_dens=df.iloc[9, 2]
spec_vol=df.iloc[10, 2]
X_dry=df.iloc[11,2]
X_wet=df.iloc[12,2]
X_eq=df.iloc[13,2]
Xs=df.iloc[14,2]
Ls= mass_MSS_start * Xs

# Initialize lists
time_sorted = [time.pop(0)]
mass_sorted = [mass.pop(0)]
time_temp = []
mass_temp = []

# Process data
for t, m in zip(time, mass):
    time_temp.append(t)
    mass_temp.append(m)

    if len(time_temp) < N:
        continue

    time_sorted.append(sum(time_temp) / N)
    mass_sorted.append(sum(mass_temp) / N)
    time_temp.clear()
    mass_temp.clear()

# Calculate X values
X = [(m - Ls) / Ls - X_eq for m in mass_sorted]

# Find the length of the longest list
max_len = max(len(X), len(mass_sorted), len(time_sorted))

# Pad the shorter lists with NaN values to make them the same length as the longest list
X += [float('nan')] * (max_len - len(X))
mass_sorted += [float('nan')] * (max_len - len(mass_sorted))
time_sorted += [float('nan')] * (max_len - len(time_sorted))

# Create the data frame
df_sorted = pd.DataFrame({'time_sorted': time_sorted, 'mass_sorted': mass_sorted, 'X': X})

# Write the dataframe to an Excel file
with pd.ExcelWriter(title+'_results.xlsx') as writer:
    df_sorted.to_excel(writer, sheet_name='Sheet1', index=False)

#GRAPH SORTED VS. UNSORTED DATA

plt.figure()
plt.plot(time, mass, ".r", label="Unsorted data")
plt.plot(time_sorted, mass_sorted, "--b", label="Sorted data")
plt.xlabel("t (h)")
plt.ylabel("m (g)")
plt.title(f"Sorted Experimental Data For ({title}) at N = {N}")
plt.legend()
plt.grid()
plt.show()
############# REGRESSION OF THE 1. DRYING PERIOD #############
x = df_sorted['time_sorted'].to_numpy().reshape(-1, 1)
y = df_sorted['X'].to_numpy().reshape(-1, 1)
r1_squared = 0
epsilon=0.996   #criteria for R_sqaured

while r1_squared <= epsilon:
    # Perform linear regression
    regression = LinearRegression()
    regression.fit(x, y)
    # Calculate coefficients and R-squared value
    k1 = regression.coef_[0][0]
    n = regression.intercept_[0]
    y_pred = regression.predict(x)
    r1_squared = r2_score(y, y_pred)
    # If R-squared value is less than "epsilon"", remove one data point
    if r1_squared <= epsilon:
        x = x[:-1]
        y = y[:-1]
# Plot data and regression line of first drying period
plt.scatter(x, y, label="Experimental data")
plt.plot(x, y_pred, color='red', label=f"Linear regression: X1(t) = {k1:.2f}t + {n:.2f}")
text_box=plt.text((x[-1]+x[0])/2,(y[-1]+y[0])/1.5,f"R_sqared={r1_squared:.3f}", fontsize=10)
text_box.set_bbox(dict(facecolor='yellow', alpha=0.5, edgecolor='gray'))
plt.xlabel('t (h))')
plt.ylabel('X1 (g/g)')
plt.legend()
plt.grid()
plt.title("Linear Regression of Experimental Data of 1. Drying Period")
plt.show()

#Convert data to list
#1+2 drying period
time_all=df_sorted['time_sorted'].tolist()
X_all=df_sorted['X'].tolist()
#1. drying period
time1 = x.flatten().tolist()
X1 = y.flatten().tolist()
X_reg1 = y_pred.flatten().tolist()
#2. drying period
time2 = [a for a in time_all if a not in time1]
X2 =[a for a in X_all if a not in X1]
#output of the critical point
crit_point=[]
crit_point.append(time1[-1])
crit_point.append(X1[-1])

############# REGRESSION OF THE 2. DRYING PERIOD #############

# Exponential function: X0 * exp(-zx)
def exponential_function(x, k2):
    return crit_point[1] * np.exp(-k2 * (x - crit_point[0]))

# Example data (replace with your own data)
x_data = np.array(time2)
y_data = np.array(X2)

# Perform curve fitting
params, _ = curve_fit(exponential_function, x_data, y_data)
k2 = params[0]

# Calculate R-squared
y_pred = exponential_function(x_data, k2)
ss_res = np.sum((y_data - y_pred) ** 2)
ss_tot = np.sum((y_data - np.mean(y_data)) ** 2)
r2_squared = 1 - (ss_res / ss_tot)

# Print results
print("k2 =", k2)
print("R-squared =", r2_squared)

# Plot the graph
plt.scatter(x_data, y_data, label='Data', color='red', marker='o')
plt.plot(x_data, y_pred, label=f"Exponential regression: X2(t) = {crit_point[1]:.2f}exp(-{k2:.2f}(t-{crit_point[0]:.2f}))", color='blue', linestyle='--')
text_box=plt.text((x_data[-1]+x_data[0])/2,(y_data[-1]+y_data[0])/1.5,f"R_sqared={r2_squared:.3f}", fontsize=10)
text_box.set_bbox(dict(facecolor='yellow', alpha=0.5, edgecolor='gray'))
plt.xlabel('t (h)')
plt.ylabel('X2 (g/g)')
plt.title('Exponential Regression of 2. Drying Period:')
plt.legend()
plt.grid()
plt.show()
############################################
X_reg2 = y_pred.flatten().tolist()
#Cosmetics of Curve 1 of the First Drying Period
factor1=2
for i in range(factor1):
    X_reg1.pop(-1)
    time1.pop(-1)
X_reg=X_reg1+X_reg2
time_fixed=time1+time2

#GRAPHING

# Individual Regression Curves
plt.figure()
plt.plot(time_all, X, "--r", label="Experimental data")
plt.plot(time1, X_reg1, "--b", label="1st drying period FIT")
plt.plot(time2, X_reg2, "--g", label="2nd drying period FIT")
plt.xlabel("t (h)")
plt.ylabel("X (g/g)")
plt.title(f"The Entire Drying Curve ({title})")

text_box = plt.text(
    (time_all[-1] + time_all[0]) / 2, (X[-1] + X[0]) / 3,
    f"X1(t) = {k1:.2f}t + {n:.2f}\nX2(t) = {crit_point[1]:.2f}exp(-{k2:.2f}(t-{crit_point[0]:.2f}))",
    fontsize=10
)
text_box.set_bbox(dict(facecolor='yellow', alpha=0.5, edgecolor='gray'))
plt.legend()
plt.grid()
plt.show()
############
#Derivative calculations#
############
dmdt_unsort=[]
dmdt_sort_1part=[]
dmdt_sort_2part=[]
#CALCULATION OF THE DERIVATIVE UNSORTED
for i in range(len(time)-1):
    dmdt_unsort.append((mass[i + 1] - mass[i]) / (time[i + 1] - time[i]))
#CALCULATION OF THE DERIVATIVE SORTED 1st PART
for i in range(len(time1)-1):
    dmdt_sort_1part.append((X_reg1[i + 1] - X_reg1[i]) / (time1[i + 1] - time1[i]))
time1.pop(0)

#CALCULATION OF THE DERIVATIVE SORTED 2nd PART
for i in range(len(time2)-1):
    dmdt_sort_2part.append((X_reg2[i + 1] - X_reg2[i]) / (time2[i + 1] - time2[i]))
time2.pop(0)

#Cosmetics of the Derivative Curve
factor2=1
for i in range(factor2):
    dmdt_sort_2part.pop(0)
    time2.pop(0)

#Curve of the Entire Derivative
dmdt_cel= dmdt_sort_1part + dmdt_sort_2part
t_cel=time1+time2

# Graph of the Derivative Curve
plt.figure()
plt.plot(t_cel, dmdt_cel, "--r", label="first derivative")
plt.xlabel("t (h)")
plt.ylabel("dm/dt (g/h)")
plt.title("First derivative of the drying curve (" + title + ")")
#plt.text(300, 0.0008, "m0="+str(parametri[0])+" [g]"+"\n"+"mf="+str(parametri[1])+" [g]"+"\n"+"W="+str(parametri[0])+" %", fontsize=10, bbox=dict(facecolor="white"))
plt.legend()
plt.grid()
plt.show()

#Drying Curve R(X)
R=[]
for i in dmdt_cel:
    R.append(-i * Ls / area)
X_reg1.pop(0)
X_reg2.pop(0)

#Cosmetics of the Drying Curve
for i in range(factor2):
    X_reg2.pop(0)
X_reg=X_reg1+X_reg2

R_drying_rate = R.copy()


Cp_pare=1.7612
Cp_vode=4.184
delta_Hizp0=2501.4
M_zraka=29 #g/mol
x=0.001 #1mm thick flat plate

delta_Hizp=(Cp_pare-Cp_vode)*T_wet+delta_Hizp0
R= -k1 * (Ls / 1000 / area) / 3600  #g/m2s
h=R*delta_Hizp/(T_dry-T_wet)
Q= delta_Hizp * R * area
ky=R/(M_zraka/1000*(X_wet-X_dry))
Dab_L=k2*x

# Plotting the drying rate graph
#Flux in m/s units
C_dry=19.635*X_dry/(1+1.611*X_dry)*1/(T_dry+273.15)
C_wet=19.635*X_wet/(1+1.611*X_wet)*1/(T_wet+273.15)
ky_ms=R*1000/18/(C_wet-C_dry)

# Plotting the drying rate graph
plt.figure()
plt.plot(X_reg, R_drying_rate, "--r", label="Drying rate")
plt.xlabel("g/g")
plt.ylabel("drying rate (kg/hm^2)")
plt.title("Drying Rate Curve (" + title + ")")
plt.legend()
plt.grid()
plt.show()

### EXPORT RESULTS ###
results=[]
results.append(R*1000)
results.append(R*1000)
results.append(Q*1000)
results.append(h*1000)
results.append(ky)
results.append(ky_ms)
results.append(Dab_L/3600)
results.append('')
results.append(crit_point[1])
results.append(crit_point[0])
results_name=['','Mass flux of 1. drying period [g/s*m2]','Heat flow of 1. drying period [W] ','Heat transfer coefficient [W/mK]', 'Mass transfer coefficient [mol/m2s]','Mass transfer coefficient [m/s]', 'Dab/L [m/s]','Critical point:','Xcrit [g/g]','tcrit [h]']

#NORMALIZED DRYING CURVE#
R_max=-Ls/area*k1
R_N=[] #Normalized drying rate
X_N=[] #Normalized moisture content
for i in R_drying_rate:
    R_N.append(i/R_max)
for i in X_reg:
    X_N.append((i-X_eq)/(crit_point[1]-X_eq))

plt.figure()
plt.plot(X_N,R_N, "--r", label="Normalized drying rate")
plt.xlabel("")
plt.ylabel("")
plt.title("Normalized Drying Rate Curve (" + title + ")")
plt.legend()
plt.grid()
plt.show()

# Define the file path
file_path = title + '_results.xlsx'
# Read the Excel file into a DataFrame
df1 = pd.read_excel(file_path)
# Add a new column with row numbers to the DataFrame
df1['Row Number'] = range(1, len(df1) + 1)
# Pad the lists with NaN values if necessary
max_len = len(df1)
results += [float('nan')] * (max_len - len(results))
results_name += [float('nan')] * (max_len - len(results_name))
# Merge the lists into a new DataFrame
df2 = pd.DataFrame({'Drying parameters results': results,
                    'Drying parameters name': results_name})
# Merge the new DataFrame with the original DataFrame
df3 = pd.merge(df1, df2, how='outer', left_on='Row Number',
               right_index=True)
# Remove the row number column
df3 = df3.drop('Row Number', axis=1)
# Save the updated DataFrame to the Excel file
df3.to_excel(file_path, index=False)
# Create a DataFrame for R_N and X_N
data = {'Normalized drying rate': R_N, 'Normalized moisture content': X_N}
df4 = pd.DataFrame(data)
# Merge the new DataFrame with the existing DataFrame
df5 = pd.concat([df3, df4], axis=1)
# Save the updated DataFrame to the Excel file
df5.to_excel(file_path, index=False)


